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The “SOS” in the title does not refer to the international distress signal, but to “solid-on-solid” 

(SOS) surface growth. The catastrophic cascades are those observed by Buldyrev et al. in interde¬ 
pendent networks, which we re-interpret as multiplex networks with agents that can only survive if 
they mutually support each other, and whose survival struggle we map onto an SOS type growth 
model. This mapping not only reveals non-trivial structures in the phase space of the model, but 
also leads to a new and extremely efficient simulation algorithm. We use this algorithm to study 
interdependent agents on duplex Erdos-Renyi (ER) networks and on lattices with dimensions 2, 3, 

4, and 5. We obtain new and surprising results in all these cases, and we correct statements in the 
literature for ER networks and for 2-d lattices. In particular, we find that d = 4 is the upper critical 
dimension, that the percolation transition is continuous for d < 4 but - at least for d 7 ^ 3 - not in 
the universality class of ordinary percolation. For ER networks we verify that the cluster statistics 
is exactly described by mean field theory, but hnd evidence that the cascade process is not. For 
d = 5 we find a first order transition as for ER networks, but we find also that small clusters have a 
nontrivial mass distribution that scales at the transition point. Finally, for d = 2 with intermediate 
range dependency links we propose a scenario different from that proposed in W. Li et al, PRL 
108 , 228702 (2012). 
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I. INTRODUCTION 

While percolation and related epidemic processes had 
appeared until recently as a mature subject that holds 
few surprises, this has changed dramatically during the 
last few years [1]. A host of new models have been 
proposed, such as percolation on growing networks [5], 
percolation on hierarchical structures [3], agglomerative 
percolation HE], explosive percolation [^, percolation 
on interdependent networks [7] , percolation where nodes 
cooperate to infect their neighbors [HHin], and percola¬ 
tion where spreading agents cooperate mm- Indeed, 
some of these models are not entirely new and are related 
to models studied since the 1970’s. The model of [3], 
e.g., can be viewed as a version of percolation on lattices 
with long range contacts mnii, while cooperative per¬ 
colation [8l-fT0] can be viewed as a variant of bootstrap 
percolation [HHni m with heterogeneous nodes m- 
But these models had previously been widely considered 
as curiosities, while only the recent developments have 
shown their wide range and wide spread applicability. 

The range of behaviors found in these models is bewil¬ 
dering. Instead of the continuous (“second order”) tran¬ 
sition with standard finite size scaling (FSS) observed in 
ordinary percolation (OP), one finds everything from infi¬ 
nite order transitions with Kosterlitz-Thouless (KT) type 
scaling [3] to first order transitions with KT type scaling 
m to first order hybrid transitions uni mini [IS], and 
- last but not least - second order transitions with com¬ 
pletely different FSS behavior [20] . 

Some of these results were obtained by simulations, 
others were derived analytically. With very few excep¬ 


tions nia, the latter rely on the fact that mean field the¬ 
ory becomes exact on random locally tree-like networks 
[2TIl23j . The latter allows very elegant treatments based 
on self-consistency equations derived using message pass¬ 
ing arguments, as e.g. demonstrated in nmm for co¬ 
operative percolation and in [24] for several interdepen¬ 
dent network models [zimmEii- But these results can 
be very deceptive - as seen by simulations mm - when 
applied uncritically to networks which are either not ran¬ 
dom or not (locally) tree-like. For such networks dynamic 
message passing [SSliSH] and the cavity method [3TJ |32] 
have been applied recently with very promising results, 
but their applicability to most of the above models is still 
far from obvious. The fact that network topology can 
change the critical behavior entirely was demonstrated 
recently m cooperative spreading agents (so-called 
syndemics or co-infections). Another warning should be 
the fact that interdependent two- and three-dimensional 
lattices do not show the first order transition found on 
random networks, but show second order transitions m, 
if all links are short range. 

All this shows that efficient methods to simulate such 
models are badly needed. While there exist very effi¬ 
cient methods for OP and for all versions of cooperative 
percolation, this is not the case for the class of models 
introduced in [7] and developed further in 
l43] . For these models which can either be viewed as 
describing cascading failures on interdependent networks 
or viable clusters on multiplex networks, the algorithms 
used so far in large simulations are extremely slow. After 
the present work was done, we learned of a recent algo¬ 
rithm [44] that is fast but highly non-trivial, and which 
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has so far not yet been applied to any real problem |45j . 
It is the purpose of the present paper to present such 
an algorithm and to use it for several different network 
topologies. As we shall see, the results are most dis¬ 
turbing, as the behavior is different for each case. This 
shows that one should be extremely cautious in apply¬ 
ing results obtained mathematically for random locally 
tree-like models to real world situations. 

In deriving the algorithm we use a mapping of the 
problem onto a solid-on-solid (SOS) [46] type growth 
model. This mapping is also of interest by itself, as it 
shows that the model has a number of non-trivial features 
that might become useful in future analytic treatments. 

In the next section we shall define the model more 
formally, discuss interpretational differences with [7] and 
the mapping onto an SOS type model, and present our 
fast algorithm. Applications to Erdos-Renyi networks 
are treated in Sec. 3A, where we shall present arguments 
that the static properties of viable clusters are indeed 
as described by mean field theory, but not the dynamics 
of cascades. Applications to regular lattices in 2 to 5 
dimensions are presented in Secs. 3B to 3E. In particular 
we shall show that the percolation transition for d = 2 
is not in the universality class of ordinary percolation, 
and that d = 4 is an upper critical dimension. Finally, 
in Sec. 3F we shall discuss the behavior on 2-dimensional 
lattices with intermediate range links. The results are 
summarized, and open problems are pointed out, in the 
final Sec. 4. 


II. THE MODEL, ITS MAPPING ONTO SOS 
TYPE GROWTH, AND THE RESULTING 
ALGORITHM 

The models studied in jT] I24U271 l32ti43] were origi¬ 
nally presented as interdependent networks showing fail¬ 
ure cascades subsequent to random removal of nodes, but 
as noted in (Ml HZ] they are more easily interpreted as 
single multiplex networks. More precisely, while in gen¬ 
eral interdependent networks some nodes depend on each 
other, these dependencies always were assumed in most 
of these papers to be pairwise and mutual. In this case 
each pair of interdependent nodes can be identified into 
a single node, and the network becomes just a duplex 
network, i.e. one set of nodes connected by two sets of 
(undirected) links. In some other papers |24l |26l [34] this 
was generalized to m networks with mutual dependen¬ 
cies among all nodes in clusters of sizes < m, where each 
cluster contains at most one node per network. In this 
case identification of all nodes in such clusters leads to 
TO-plex networks. In the following we shall only consider 
the case m = 2, i.e. duplex networks or pairwise mutual 
dependencies. 

Another slight but important shift was made by [24l 
EZ] when they noted that the model can be formulated 
as a self-consistency problem without any reference to 
node removals (“damage”) and to cascades triggered by 


them. Rather, as also pointed out in [JB], it describes 
(for m = 2) the case where each node needs two essential 
supplies for being active. More precisely, we assume that 
there exist a source node that supplies both resources, 
and that the resources can be transported only on active 
nodes. It is this latter interpretation which we adopt: 
We usually consider a cluster C of nodes as viable, if any 
two of its nodes are connected by paths on both sets 
of links - where we also demand that both paths are 
entirely confined to C. In that case, any node in C can 
be supplied with both resources, if the source node is 
also in C. Conversely, if any node in C can obtain both 
resources, i.e. if any node in C is doubly connected to 
the source, then (by the assumed undirectedness of the 
links) also two arbitrary nodes in C are doubly connected. 
Notice that in this way we do not consider the networks 
themselves as failed or intact, but we consider only the 
activities on the networks as present or absent. To use the 
main example used in |7]: If there is a power failure, the 
electricity network itself might still be perfectly intact, 
it is just the activities on the power grid and on the 
associated computer network that have gone down. 

The next difference with the bulk literature on inter¬ 
dependent networks is that we do not in general delete 
nodes, but we consider the fate of viable clusters as we 
change the densities of links. On ER lattices, decimating 
nodes is just equivalent to renormalizing the connectivity 
and thus equivalent to decimating links [27]. On regu¬ 
lar lattices this is not true, but keeping all nodes just 
reduces the number of control parameters by one. For 
two-dimensional (2-d) lattices we checked, however, that 
decimating both sites and links leads to the same univer¬ 
sality class. 

A last difference with [7] is that we consider not only 
the largest viable cluster on the entire network, but we 
study all viable clusters, in agreement with [33]. As we 
shall see, this comes with no additional effort. Indeed, 
only by studying all viable clusters we can verify with 
our algorithm which one is the largest one. Again we 
consider this as more realistic. To stay with the electric 
breakdown example: Assume we have a power failure in 
central Italy which disconnects the north from the south. 
This would be considered by [Z] as catastrophic, since 
there would no longer exist a giant viable cluster. But 
as long as there local viable clusters around Milano and 
Napoli, people there would be perfectly happy. 

In spite of all the differences with [7] mentioned in this 
section, notice that we are mathematically still dealing 
with precisely the same class of models, and expect the 
same type(s) of phase transitions. 

Our algorithm works for any type of duplex network 
with undirected links, although the mapping onto an SOS 
growth model requires of course strictly spoken that we 
deal with a planar lattice (since SOS models are mod¬ 
els for 2-d surfaces). We shall use therefore a language 
adequate for this special case, although it should be un¬ 
derstood that everything in the following applies also the 
general case. We thus have a set of N points which 
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are partially connected by two sets of bonds (“red” and 
“green”, say) between neighboring sites. Bonds of either 
color are placed randomly and independently, with prob¬ 
ability q G [0,1]. 

Our algorithm has two ingredients. In the first part we 
pick a seed (or “source”) site and find the largest viable 
cluster connected to this seed. In the second part we 
repeat this for all possible seeds. 

(a) Finding the largest viable cluster C attached 
to point V. This is simply done by alternatively per¬ 
forming “epidemic” or “Leath-type” m spreading pro¬ 
cesses on the red and green bonds, each starting from 
site i (this can be done breadth or depth first; we actu¬ 
ally used breadth first. We also assume for definiteness 
that we start with the red bonds). We do not fix bond 
occupancies ” on the run”, but we rather determine them 
before we start with the first epidemic. We follow the 
spreading until it dies due to the finiteness of the lattice, 
which gives us a first cluster Ci. Since Ci is connected 
to the seed site but not necessarily doubly connected, 
we know that the largest viable cluster attached to site 
i must be contained in it, C C Ci. Therefore, when we 
generate the second epidemic using the green bonds, we 
restrict ourselves to Ci, generating thereby a cluster 
with C C C2 C Cl . As we proceed alternatingly, we gen¬ 
erate thus a chain of nested sets 

C... CC/, CC,,_i... CCi. ( 1 ) 

Since the lattice is finite, this must stop at a finite 
value of h which we call hi. It is easy to see that Chi 
is equal to C. If not, then Chi would contain at least 
one point which is not connected to the seed by paths 
of either color, but in that case Chi+i would be strictly 
larger than C/i^, which is in conflict to our assumption 
that the iteration stops at hi. 

In the following, we shall call each epidemic starting 
from i a “wave”, and h its height. Since the supports Ch 
of successive waves are nested, we obtain in this way a 
landscape with a single mountain of peak height hi, and 
“terraces” of heights h < hi (see Fig. I). Each terrace is 
just 

ACh=Ch\Ch+i. ( 2 ) 

The part of the lattice which was not even touched by the 
first wave is called the zeroth terrace with height h = 0 

m- 

(b) Finding viable clusters attached to other 
points: After we have obtained the largest cluster at¬ 
tached to i we could pick another site j (either randomly 
or by going systematically through the lattice), and re¬ 
peat the same process. To find the largest viable cluster 
on the entire lattice Cmax) we would then have to repeat 
this for all sites. But this would be extremely time con¬ 
suming. We are interested in the phase transition where 
Cmax becomes macroscopic in the limit N ^ 00 . This 
happens at a critical value of q which is far above the 
critical value for single epidemics. Thus for all seeds the 



FIG. 1. (color online) panel (a): Clusters Ci (all non-white 
points) to C5 (dark blue) attached to the site i indicated by 
the black dot. Each color corresponds to one terrace, with 
light blue < green < orange < magenta < red indicating in¬ 
creasing heights; panel (b): six large clusters of height 2, all 
attached to points in the cluster Ci shown in panel (a). Notice 
that all these clusters are fully contained in Ci, and no two 
clusters overlap partially. They are either disjoint, or one is 
a subset of the other. Notice also that some of these clusters 
touch each other, indicating the possibility that some adja¬ 
cent terraces may have equal heights. The control parameter 
is set such that ordinary percolation would be supercritical, 
but mutually interdependent percolation of viable clusters is 
subcritical. 


first few waves will most likely be huge, and most of the 
CPU time is spent for finding out clusters which overlap 
so much that they are nearly identical and cover nearly 
the entire lattice. The resulting CPU time would then 
be roughly O(Af^). 

Fortunately, there is no need to use this brute force 
algorithm. Let us consider the system of waves and ter¬ 
races for the second seed j. Since the entire lattice is 
covered by terraces from the first seed i, we can assume 
that j is on a terrace of some height h G {0,1 ,...hi}. 
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If the height were equal to hi, j would be part of the 
first cluster. Thus we can assume that h < hi. If /i = 0, 
then j is not connected to i by the red bonds, and thus 
any viable cluster attached to j must be contained in the 
zeroth terrace ACq. Otherwise {if 0 < h < hi) we can 
use 

Lemma 1: If h > 0, the first h waves starting from j 
cover exactly the same clusters Ch' (with h' < h) as the 
waves starting from i. 

Proof: The proof is by induction and uses the fact 
that each wave just covers all sites on a given network 
that can be reached from the seed, where each of these 
networks is the subgraph of the original network reached 
by the previous wave. First of all, it is clear that the 
lemma holds for h = 1, because Ci is the set of all points 
on the original network that are connected to the seed 
by red bonds. Since the two seeds are connected by red 
bonds due to the assumption that h > 0, any point k 
connected to i must also be connected to j and vice versa. 
Let us now assume that for some h' < h the clusters 
attached to i and j are the same. Then we can use exactly 
the same argument, just replacing the original network 
by Ch' and using the appropriate color of the bonds. 

Thus we do not need to follow the first h waves starting 
from j, and we can immediately start with wave h + 1. 
In doing so, we can use also another important simplifi¬ 
cation because of 

Lemma 2: Assume that point j is on a terrace of 
height h < hi. Then the entire next wave is confined to 
this terrace. 

Proof: For the proof we have to show that the next 
wave neither spills over to the lower terrace, nor to the 
higher. 

For the hrst part we can assume h > 0. Assume now 
that the wave actually does spill over, i.e. there exists a 
point k which is connected to j but is on a lower terrace. 
This point would not be connected to i on any terrace 
with h' >= h, while j is. This cannot be. 

For the second part we assume for the moment that the 
wave spills over to a higher terrace, i.e. there exists a 
point k which is connected to both k and i. But since j 
is not on this higher terrace (i.e. is not connected to i), 
this cannot be either. 

Thus when building the terraces for site j we can re¬ 
strict ourselves to waves for which all boundaries of the 
terraces generated by the first seed form obstacles which 
cannot be crossed. It is easy to see that this generalizes 
also to all later seeds, when the landscape is made up be 
any number of terraces and boundaries between them: 

Proposition 1 : If a seed j happens to be on any ter¬ 
race with a height h which is not a local maximum (i.e., 
j is not in a locally maximal viable cluster), then all pre¬ 
vious wave boundaries form obstacles for all avalanches 
starting at j which have to be followed in order to get the 
viable cluster attached to j. 

The proof uses the same arguments as above. We just 
have to realize that successive waves that started from 
later and later seeds are just distinguished by more and 


more restricted subsets of the original network to which 
they are confined. Whatever these subsets are, the above 
proofs go through without modihcation. 

It should now be clear why we call this procedure an 
SOS type growth process. As in any SOS model, we have 
a rough landscape and the growth of this landscape pro¬ 
ceeds by localized events, each of which builds a hierarchy 
of nested terraces. As in surface diffusion problems with 
strong Schwoebel barriers, these events cannot spill over 
the boundaries set by previous events. 

Notice that neighboring terraces belonging to different 
seeds can have the same height. Also in this case, their 
boundaries cannot be crossed by later waves from other 
seeds. To implement this restriction we either use an 
additional marker for each site which tells us the seed 
to which it belongs; or, alternatively, when each wave 
is finished, we cut all bonds connecting the wave with 
its complement. Both methods allow us to prevent the 
waves from spilling over the terrace boundary, without 
modifying the dynamics inside any terrace. Also, when 
we start with a new seed j, we choose the color of the 
first wave according to the landscape height h at this 
seed, and we count subsequent heights by adding to this 
h. Since now every site i is “infected” precisely by hi 
waves, it follows that the algorithm has time complexity 
N{h), where {h) is the average height. This estimate 
holds for arbitrary graphs, provided all node degrees are 
bounded. 

Before we show numerical results, we present also a 
second proposition which we did not find useful for nu¬ 
merics, but which could be very useful for mathematical 
treatments. 

Proposition 2: The final landscape and the system of 
all terrace boundaries is independent of the sequence by 
which the seed points are chosen. 

Thus the landscape is not a property of the realiza¬ 
tion of the algorithm (which involves an arbitrary choice 
of going through all seeds), but is an inherent property 
of the network. In essence, it says that the growth of 
the surface is Abelian in a similar way as the Bak-Tang- 
Wiesenfeld sand pile model is Abelian m- 

For the proof we can use the fact that any permutation 
of sites can be written as a product of pairwise transposi¬ 
tions. We thus need to show only that it does not matter 
whether we first take seed i and then seed j or the in¬ 
verse. But this follows from the arguments in the proofs 
of Lemma 1 and 2. 

We should add that we checked both propositions also 
numerically, finding perfect agreement. 

Finally, we should point out how our sequences of 
waves are related to the failure cascades of [7] and to 
other algorithms proposed for them. One essential dif¬ 
ference is that we do not demand that the cluster sur¬ 
viving the cascade is always the maximal one. Rather, 
the maximal cluster is determined at the end, when all 
viable clusters are known. This simplifies the algorithm 
of course considerably. Also the complexity of the algo¬ 
rithm of [U] is related to the fact that they always follow 
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the largest viable cluster. This precaution is not taken in 
the algorithm of [63], where it is assumed that a cluster 
which starts large at the beginning of the cascade is not 
overtaken later by one which starts smaller. This is true 
for sparse tree-like graphs (whence their algorithm gives 
correct results), but it is not true in general. Related to 
this is the absence of the notion of ‘seed’ or ‘source’ in 
these algorithms. While in our algorithm each ‘cascade’ 
gives just one cluster arbitrarily picked by the seed, the 
cascades studied in |7] are constructed so as to lead al¬ 
ways to the largest cluster. As a consequence, the average 
cascade length in |7] is in our interpretation essentially 
the average height of the largest viable cluster, which is 
in general (but not always?) an upper bound for the 
average surface height. 

A last difference between our algorithm and that of 
[3 is that we always start our ‘cascades’ from the full 
(undamaged) system, while the cascades of |7] are sup¬ 
posed to be triggered by (small and successive?) dam¬ 
ages, starting from an already partially damaged system. 
On the other hand, our algorithm shares with that of |7] 
that it does not deal - in contrast to what is suggested in 
[7] - with any real dynamics. Both algorithms deal with 
pseudodynamics in a fictitious time, and should not be 
confused with actual cascading processes going on in real 
time. 


III. APPLICATIONS 
A. Erdos-Renyi Networks 

We first apply our algorithm to duplex ER networks 
where both layers have the same average degree z = (fc). 
In this case it is known from HlllTI that the generating 
function formalism of mmn gives the exact threshold 
and the exact dependency of the order parameter on (fc). 
Thus it can be used to test the accuracy of the algorithm. 
On the other hand, we shall also use it to find the aver¬ 
age height of the landscape and the average height of its 
maxima (i.e., of the largest clusters). The latter is just 
the average cascade life time in the interpretation of [7]. 
Since during the cascades the clusters are supercritical, 
it is not a priori clear whether the mean field theory ar¬ 
guments used in |7] to derive a scaling law for this life 
time are exact. To make a short summary, we will find 
that the algorithm works perfectly, but the mean field 
arguments of |7] for the life time seem to be to be only 
approximate. This implies that also the arguments of 
[7] in favor of the percolation threshold value have to be 
taken with care - although they give the correct result 
- because they rely on the assumed mean field cascade 
dynamics. 

We consider networks with iV = 2"* nodes, with m = 
4,... 26. The average degree was varied in the range 
1.9 to 3.2. The critical point is known in this case to 



FIG. 2. (color online) Order parameter (density of giant vi¬ 
able clusters) for ER networks plotted against the average 
degree 2 . The continuous curve is the prediction of Eq. §, 
while the points are from simulations with N ranging from 
2^"^ to 2^®. Eor S < 0.511700... the theoretical curve is strictly 
vertical, while it has square-root behavior above. 



FIG. 3. (color online) Differences between the data points 
and the theoretical curve in Fig. 2. To enhance the region 
z ~ Zc, the data are plotted against z — Zc on a, logarithmic 
axis. 


be at Zc = 2.45540748.... The number of realizations 
simulated at Zc varied from « 10® for N = 1024 to > 600 
for N = 2^^. 

The density of the giant viable cluster on an infinite ER 
network is given by the largest solution of the equation 

m 

F{S) = S-{l-e-^^)‘^ =0. (3) 

It vanishes for z < Zc, and has both a jump and a square 
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FIG. 4. (color online) Differences between the data points for 
z = Zc and the theoretical value 0.511700 towards which they 
should converge for N ^ oo. We see indeed a straight line on 
a log-log plot, which leads to Eq. ([^. 

root singularity at z = Zc, 

S = 0.511700 -I- aVz — Zc + ... for z > Zc- (4) 

For finite N there are of course corrections, the analytic 
form of which is not known. 

Numerical results for Sn, the average relative size of 
the giant cluster on a graph of size N, are given in Figs. 2 
to 4. Actually, in these plots we show (in contrast to 
analogous plots in the next subsections) Sjy conditioned 
on realizations which do have a giant cluster. Notice 
that all other clusters are extremely small (for N > 
we found only clusters of sizes 1 or 2 apart from the giant 
one, and even for N = 1024 the largest non-giant clusters 
had sizes < 8), thus it is straightforward to identify giant 
clusters, except for N < 10^. In Fig. 2 we show Sn 
versus z. Since this gives a perfect agreement with Eq. @ 
except for z extremely close to Zc, we plot in Figs. 3 and 4 
differences between analytical and numerical results for 
a more sensitive comparison. From Fig. 3 we see that 
the finite-iV corrections are indeed not monotonic, as one 
might have suggested from Fig. 2, but they decrease fast 
with N. For N = 2^®, the deviations from the theory 
for iV = 00 are < 10“"*^ for all z except extremely close 
to Z 4 , which is just the level of statistical errors. Thus 
we verified with very high precision that our algorithm 
gives results in agreement with Eq. ([^ not only near the 
critical point, but also for all z > Zc- The agreement at 
z = Zc is checked in Fig. 4, where the differences ASn = 
Sn — Soo are plotted against N. We see a perfect scaling 
law 

ASn with a = 0.23(1). (5) 

To our knowledge there is no theoretical prediction for a. 
For z > Zc, Fig. 3 is consistent with ASn ~ N~°‘(l)[{z — 


Zc)A^^/^], but the data are too noisy make a strong claim 
for it. 



(z - Zc ) “ 


FIG. 5. (color online) Pn{z), the probability to have a giant 
viable cluster on an ER network of size N, plotted against 

(z-Zc)Nh2. 

There also seems to exist no prediction for the behav¬ 
ior of Pn{z), the chance that there exists a viable giant 
cluster on a network of size N. In analogy to ordinary 
percolation we expect that this is a step function in the 
limit N = oo, i.e. Poo{z) = — Zc), but the behav¬ 

ior for finite N is non-trivial. In Fig. 5 we plot Pn{z) 
against (z — Zc)IV^/^. We see that all curves become es¬ 
sentially parallel, i.e. we have a data collapse apart from 
a correction to scaling that leads to a shift of the effective 
transition point towards higher z as N increases. 



FIG. 6. (color online) Average “surface” heights (i.e. cascade 
life times) for three network sizes, plotted against z. While the 
collapse of the three curves for z <C Zc is easy to understand, 
the behavior at larger z is more interesting. 

“Surface” heights (i.e. average cascade life times) mea- 
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sured in these simulations are discussed in the next fig¬ 
ures. In Fig. 6 we show for three different network sizes 
how the average surface heights, averaged over all nodes, 
depend on z. As expected, they are small for z Zc 
and become quickly independent of N : For small z there 
are only small clusters to start with, and the cascades 
terminate quickly. For z ^ Zc there are only small holes 
in the viable clusters and they also are qnickly identified, 
whence the cascades terminate soon also in this regime. 
It is less obvious why the heights do not seem to be¬ 
come independent of as z becomes very large. The 
most conspicuous feature of Fig. 6 is, however, the sharp 
peaks developing at z = Zc as iV increases. Figures 7 to 
9 deal with several aspects of these peaks. 

Heights of the peaks (more precisely the average sur¬ 
faces heights at z = Zc; the peaks occur slightly to the 
left of Zc) are shown in Fig. 7. In the same figure we 
also show, in addition to the average surface heights, the 
(averages of the) surface peaks, i.e. the highest points in 
the “landscapes”. At least for ER networks (where we 
have no large viable clusters apart from the giant ones) 
we expect these peak heights to be precisely the average 
durations of the cascades studied in [7] (see also @3]). We 
see from Fig. 7 that the average heights show a perfect 
scaling. 


{h) ^ , 7 = 0.280(1). ( 6 ) 

In contrast, the peak heights in the landscape seem to 
increase with the same asymptotic power law, bnt with 
large finite-A^ corrections. The latter are presumably the 
usual bias corrections in extremal properties evaluated 
over finite domains. One has similar (logarithmic?) cor¬ 
rections e.g. in the average size of the largest cluster in 
subcritical (ordinary) percolation on finite lattices. 

Eq. Q is in clear contradiction to the result /ipeak 
obtained in [7] by mean field arguments which do 
give also the correct Eq. (|^ for the cluster size and the 
exact value of Zc- This is rather puzzling. One might try 
to explain it by noting that the value of Zc depends only 
on clusters which are barely viable, and thus the local 
tree-likeness of sparse ER networks should be sufficient - 
while the cascade dynamics deals always with supercrit¬ 
ical clusters. But then it is not clear why Eq. (|^ should 
be correct also for z > Zc, which we took great pain to 
verify. Notice that the violation of the mean field pre¬ 
diction seen in Eig. 6 is not related to the observations 
in |43) . In that paper it was shown that an apparent vi¬ 
olation of mean field theory is observed, if one does not 
use the true critical value of Pc in the analysis, but an 
effective one that changes from realization to realization. 
In our analysis, Pc is always the true critical point. In 
any case, our value of 7 disagrees both with the mean 
field exponent 1 /4 obtained in [7] and with the exponent 
1/3 proposed in |43| . which was based on simulations of 
much smaller systems with much lower statistics. 

The “wings” of the peaks seen in Eig. 6 are studied 
more closely in Fig. 8 , where we plotted the heights (on a 
log scale) against log |z — Zc|. For the subcritical regime 
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FIG. 7. (color online) (a) Average “surface” heights at z = 
Zc plotted against networks network sizes on a log-log plot, 
and average values of the peaks in each surface. Both curves 
clearly disagree with the scaling oc proposed in [7], but 
give rather an exponent 0.280(1); (b) The same data, but 
divided by 


z < Zc we clearly see a power law with exponent 1 / 2 , 
in perfect agreement with [43] (notice, however, that we 
could not confirm the other scaling laws suggested in 
|43) 1. We did not try to derive this analytically, but do¬ 
ing this should not be difficult. For z > Zc the situation 
is less clear, but it seems that for large N a power law 
with exponent 1/3 emerges, even if corrections are large 
for finite N. 

The two panels of Fig. 9 show finally two attempts to 
produce a data collapse for the entire critical region. In 
panel (a) we plotted {h)/N'^ against (z — Zc)N‘^^. This 
gives us a perfect data collapse subcritically. In the crit¬ 
ical and supercritical regions the plot does not look too 
bad either, but a closer inspection shows that it is in¬ 
deed unacceptable. In contrast, in panel (b) we plotted 
the data against (z — Zc)N^'^ ■ This time the collapse is 
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FIG. 8. (color online) Log-log plots of (h) against \z — Zc\- 
The straight lines indicate power laws which we claim to hold 
in the limit N ^ oo. 

very poor in the wings, but we get a very good collapse 
in the critical region - apart from the same shift of the 
effective critical point also seen in Fig. 5. Notice that the 
variables on the x-axis in Figs. 5 and 9 are the same. 


B. 2-dimensional Lattices 

We next study 2-d square lattices of size L x L, with 
32 < L < 32768 and with helical boundary conditions 
[M] . The two sets of links are different (randomly chosen) 
subsets of nearest neighbor bonds. For each of them we 
include a bond with probability p and exclude it with 
probability 1 —p. In addition, we can also make a fraction 
q of sites inaccessible. For p = 1 and 0 < q < 1 the 
red and green links coincide, and we simply deal with 
site percolation with control parameter q. On the other 
hand, if p < 1 and q = 1, the sets of red and green 
links overlap. Each one of them would define a bond 
percolation problem, and the viable clusters are subsets 
of intersections of red and green percolation clusters. 

For p < 1 and q < 1, the single color problems would 
correspond to mixed site-bond percolation. We studied 
also that case, because the problem was originally formu¬ 
lated in [7] as a robustness problem for network break¬ 
down under successive node removal. In that interpre¬ 
tation one needs p < 1 for having a multiplex network, 
and (? < 1 due to node removal. Actually, however, the 
problem just becomes more complicated by having both 
p and q, without adding any conceptual advantage. If 
one wants to interpret results obtained with q = 1 (i.e. 
with an undiluted lattice) in terms of robustness against 
damage, then one simply has to implement the damage 
not as site removal, but as bond removal. 

In the following we shall present results for q = 1 and 
varying p (because this is conceptually the simplest), and 
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FIG. 9. (color online) Data collapse plots of re-scaled surface 
heights {h)/N^ against {z — Zc)N^'' with 2'y = 0.56 (panel 
a) and against {z — Zc)N^'^ (panel b). Neither plot is fully 
satisfactory, but in spite of the bigger overall discrepancies we 
claim that panel (b) gives the correct collapse in the critical 
region. 


for p = 0.6 and varying q. The latter was done in order 
to compare with the results of [za EH [52]. It leads to 
Qc ~ 0.96. In neither case the exact critical point is 
known. The total number of realizations were in both 
cases between > 10® for the smallest, and several thou¬ 
sand for the largest lattices. 

Results for p = 0.6 are shown in Figs. 10 to 12. In 
Fig. 10 we show a data collapse similar to the corre¬ 
sponding plots in EH EH- We see a perfect collapse, 
if we use qc = 0.960512(4) and ix = 1.2. These values 
agree roughly with those proposed in but are very 
different from those in m- In particular, we would ob¬ 
tain a very bad collapse (for any value of ix), if we would 
use qc = 0.9609 as proposed in m- Also, we would ob¬ 
tain a very bad collapse (for any value of qc ), if we would 
take for v the value 4/3 as in ordinary percolation. We 
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FIG. 10. (color online). Data collapse plot for the size of the 
largest viable cluster on 2-d lattices with p = 0.6. The only 
perceptible deviation from a perfect collapse is at very large 
values of q, where the data for L = 256 and for L = 512 are 
systematically too high, due to finite-size corrections. 


FIG. 12. (color online) Log-log plot of against L, 

for eight values of q close to qc- The power 91/48 = 1.8958 ... 
was chosen, because it is the fractal dimension of the incipient 
giant cluster in critical OP. Thus any deviation of the critical 
curve from a horizontal line indicates that the present model 
and OP are in different universality classes. 



FIG. 11. (color online). Log-log plot of S against q — q^, 
demonstrating the power law S' ~ (5 — qc)^■ 


should point out that the largest lattice sizes studied in 
[5T] were only 3000 x 3000. 

Supercritical data are plotted as S against q — qc on 
a log-log plot. For large L and q not too close to qc we 
expect a power law 


Finally, we show in Fig. 12 a log-log plot of S against 
L. To make the plot more signihcant, we actually divides 
S by the power of L expected for OP, so that the critical 
curve would be horizontal if the present model were in 
the OP universality class, as claimed in m- The dashed 
straight line indicates the power law consistent with the 
previous two hgures. 


(n 

CO 

> 



S'^iq-qcf. (7) 

This is indeed nicely confirmed, with /3 = 0.163(2). No¬ 
tice that this is perfectly consistent with the scaling re¬ 
lation D = 2 — /3/jz, but is again different from the value 
for OP. Notice also that the power law holds only for 
q — qc < 0.01. For larger distances from the critical point, 
corrections would be needed. 


FIG. 13. (color online). Data collapse plot for the variance 
of the size of the largest viable cluster on 2 -d lattices with 
q = 1 . 0 , plotted against rescaled values of p — Pc- 

Essentially the same results were obtained when we 
kept q = 1 fixed (i.e., no sites were deleted), and p was 
allowed to vary. This time we studied only lattices with 
sizes up to 16384 x 16384, but with higher statistics. The 
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critical point is now pc = 0.576132(5). All critical expo¬ 
nents are compatible with these given above, and have 
roughly the same error bars. The main difference with 
the case p = 0.6 is that now corrections in the far super¬ 
critical region do not make the curves in the log-log plot 
of S against p — pc turn up (as in Fig. 11), but make then 
turn slightly down. To avoid duplication we do not show 
the dame plots as for p = 0.6, but we do show a data 
collapse plot for the variance (“susceptibility”) of S. 



FIG. 14. (color online) Average SOS landscape heights for 
2 -d lattices with q = 1, plotted against p. 

Average landscape heights are shown in Fig. 14. We 
see a behavior rather different from that seen in Fig. 6 
for ER graphs. While the latter showed a very sharp 
peak at the critical point, whose height and sharpness 
both increased with powers of the system size, we have 
now a very wide bump far in the supercritical region, 
whose width does not seem to depends on L at all, and 
whose height increases logarithmically with L. The most 
conspicuous feature of Fig. 14 is that the slopes at p = Pc 
become increasingly steeper with L, roughly as 

^(/i)~L°'^ at p = Pc- (8) 

dp 

In contrast to the ER case, where we could find a de¬ 
cent scaling behavior of {h) near the critical point, we 
were unable to find any convincing scaling in the present 
case. 

Before leaving this subsection, we should mention that 
all scaling laws and critical exponents found in this sub¬ 
section were also found in a very different microscopic 
realization of the 2-d lattice model discussed in subsec¬ 
tion III.F. 

C. 3-dimensional Lattices 

For lattices with d > 2 we made only simulations with¬ 
out site decimations, i.e. we always used q = 1. 


In HZl, simulations of very small lattices had suggested 
that the transition is also for 3-dimensional lattices con¬ 
tinuous and not in the OP universality class. This is 
important, since it could be argued that the transition is 
continuous for d = 2 because there the two sets of links 
necessarily overlap strongly, and in that case the tran¬ 
sition can be continuous even on ER networks [25) . In 
three dimensions, this overlap is much reduced. 



L 


FIG. 15. (color online) Log-log plot of Li^S /against 
L, for eight values of p close to Pc- Here, Dqp = 2.523(1) 
is the fractal dimension in three-dimensional OP |53l I54| . 
Indeed, since we used helical b.c. where the number 
N of lattice sites was always a power of 2, some of 
the lattices were not strictly cubes. In these cases, L 
was nevertheless defined as L = , even if this was 

not an integer. The values of p are (from top to bottom) 
0.357, 0.3568,0.35673,0.35670,0.35668,0.35665,0.3566, 0.3565. 

Our present simulations show that the transition is in¬ 
deed continuous, but the deviations from OP scaling are 
much smaller than found in m- In Fig. 15 we show 
a log-log plot of rescaled infinite cluster sizes against L 
for eight values of p close to Pc- Here L is defined as 
L = where N is the number of lattice sites. Since 

we used helical b.c. where N was always a power of 2, 
the so defined L is not always an integer, and the lattices 
are not always exact cubes. This has however no notice¬ 
able effect on the scaling. Since we want to compare with 
OP, where L^S ^ with Dqp = 2.523(1) [55115^ . 
we multiply the data with the corresponding power of L. 
Thus we expect one curve (the critical one) to be horizon¬ 
tal, if and only if the present model is in the universality 
class of OP. We see that the most straight curve has a 
slightly negative slope. This gives Pc = 0.356707(4) and 

Df = 2.475(8) (9) 

(the estimate of [27] was 2.40(1)), which would nominally 
indicate that the present model is in a different universal¬ 
ity class. As we shall see, however, there are important 
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corrections to scaling. Thus the error in this estimate 
may be severely underestimated. 



FIG. 16. (color online). Data collapse plot for the size of the 
largest viable cluster on 3-d lattices. In contrast to the 2-d 
case we now have huge corrections to scaling in the supercrit¬ 
ical region. 
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FIG. 17. (color online) Average SOS landscape heights for 
3-d lattices, plotted against p. 

These corrections to scaling are most clearly seen in a 
data collapse plot, see Fig. 16. In contrast to the analo¬ 
gous plot for d = 2 in Fig. 10, we see now huge deviations 
in the supercritical region. Nevertheless, in the critical 
and under critical regions the collapse seems to be very 
good. Notice that we used here the values for Dj and 
Pc obtained from the previous plot. For the correlation 
length exponent we obtain v = 0.80(3), which is signifi¬ 
cantly lower that the value 0.8734 obtained fro OP [54]. 
Hyper scaling finally gives /3 = (d — Dj)!^ « 0.42, in 
agreement with OP. In summary it seems that correc¬ 
tions to scaling were severely misjudged in HZ], and we 


cannot exclude that the model is in the OP universality 
class for d = 3. 

A reason for it to be not in the OP universality class 
is, of course, the divergence of the surface height when 
L —>■ oo. This is clearly seen from Fig. 17, where we 
find the same logarithmic increase as in d = 2, and a 
power law for the slope at p = Pc as in Eq. <§ , but with 
exponent « 0.82. 

D. 4-dimensional Lattices 



FIG. 18. (color online) Average SOS landscape heights for 
4-d lattices, plotted against p. 

As a first result we show the average surface heights, 
see Fig. 18. Their behavior is very similar to those in two 
and three dimensions: For any N = their maxima oc¬ 
cur in the supercritical region, and their heights increase 
logarithmically with L. Also, the slopes d{h)/dp diverge 
for L —^ oo with a power law, again with power slightly 
smaller than one (our best estimate is « 0.9). 

In contrast to this, the behavior of the order param¬ 
eter is completely different. In Fig. 19 we show plots 
of S versus p, for lattice sizes with N = 2^®, 2^®,... 2^®, 
corresponding to L = 16,... 128. We see clearly a con¬ 
tinuous transition, but with order parameter equal to (or 
at least very close to) (3 = 1. This would indicate that 
d = 4 is the upper critical dimension, in striking con¬ 
trast to OP, where d„ = 6. An attempt to collapse the 
data according to standard FSS leads to Fig. 20, and to 
Df « 8/3 and « 3/4. These would not correspond to 
any known mean field theory, but we should be aware 
that we should expect logarithmic corrections, if indeed 
du = 4. The poor quality of the collapse might indi¬ 
cate that such corrections are indeed present, but with 
our present data we have no chance to estimate them in 
detail. 

Finally, we show in Fig. 21 the integrated mass distri¬ 
bution on a large lattice for p « pc- If usual FSS would 
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FIG. 19. (color online) Average density of the largest viable 
cluster on 4-d lattices with N sites, plotted against p. The 
sharpness of the kink at p « 0.2715 increases with N. 


FIG. 21. (color online) Log-log plot of the integrated mass 
distribution of 4-d viable clusters at L — 128 and p = 0.27145. 
The straight line has slope —3/2 as predicted from standard 
scaling laws, accepting the exponents found in Fig. 20. 
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FIG. 20. (color online) The same data as in Fig. 19, but 
plotted such that the collapse onto a single curve near the 
critical point. 


E. 5-dimensional Lattices 



hold, we would expect a power law 

p{m) ~ m^~'^ (10) 

with 


FIG. 22. (color online) Average density S of the largest viable 
cluster on 5-d lattices with N sites, plotted against p. Notice 
that we averaged here over all runs (even if they did not have 
a giant viable cluster), while the averages in subsection III.A 
(e.g. Fig. 2) were done only over runs that did have a giant 
cluster. 


= 1 + d/Df . 


( 11 ) 


Inserting here our numerical result Df = 8/3 gives 
T = 2.5. The data show a very rough power law, with 
large deviations at intermediate mass values, but the dis¬ 
tribution agrees surprisingly well with Eq. (101 for large 


masses. 


If the present problem represents indeed a new univer¬ 
sality class with upper critical dimension = 4, then 
the behavior in d = 5 should be essentially the same as 
for d —>■ oo, and in particular as on ER graphs. To ver¬ 
ify this we show first (in Fig. 22) the order parameter S 
(the density of the largest viable cluster) as a function 
of p. We see that the curves for different sizes N inter¬ 
sect indeed in a single point, with pc = 0.22614(1) and 
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Sc = 0.31(1). This jump of S' is a clear sign of a dis¬ 
continuous transition, although the rise of S for p > pc 
indicates that the transition is also, as for ER networks, 
hybrid. 



P 


FIG. 23. (color online) Average SOS landscape heights for 
5-d lattices, plotted against p. 

The same conclusion is reached by looking at surface 
heights (Fig. 23), which display the same sharp peak as 
for ER networks. Also, the heights of the peaks increase 
with roughly the same power of the system size (/imax ~ 
7V0-30(2)^^ not logarithmically as on lattices with d < 4. 



FIG. 24. (color online) Log-log plot of the integrated mass 
distribution of finite 5-d viable clusters at A = 2^® (L ~ 37) 
and p = 0.22615. The straight line indicates the scaling of 
finite but large viable clusters. 

On the other hand, the statistics of finite viable clus¬ 
ters is markedly different from that for ER networks. 
There, apart from the giant viable cluster all other clus¬ 
ters are of size one in the TV —> oo limit [7], and for 


the system sizes studied in the present work most of 
them have sizes one or two. On five-dimensional lattices, 
however, the distribution of finite viable clusters is non¬ 
trivial. The average mass of the second-largest clusters 
peaks for Pc — P ^ and the peak height increases 

roughly as (both of these exponents are very crude 
estimates). Moreover, the mass distribution of the finite 
clusters shows rough scaling, p{m) ~ dX p = Pc 

(see Fig. 24). 


F. Cross-over to Mean Field Theory via Long 
Range Links on 2-d Lattices 

1. General Remarks: Long Range Connectivity vs. Long 
Range Dependency 

In general, there are two main ways how the anomalous 
behavior of critical phenomena seen on low dimensional 
lattices can cross over to mean field theory. One is by 
increasing the dimension, the second is by making the 
theory more and more non-local. In the present model we 
found that the transition via dimension increase is very 
non-standard, by showing first a transition to a mean 
field critical point (at c? = 4), which is then replaced by 
a first order transition at d > 4. In view of this, it is 
of interest to check what mechanism(s) are at play when 
one makes the model more and more non-local. 

In the interpretation as multiplex networks (where any 
interdependent node pairs are joined into a single node) 
adopted in this paper, the most natural way to do this 
is by increasing the lengths of the (“connectivity”) links. 
In analogy with standard critical phenomena one would 
expect that the model stays in the same universality class 
as long as the characteristic length of the bonds are fi¬ 
nite, and crosses over to mean field theory only when this 
length becomes infinite, e.g. by having power-behaved 
link length distributions (see, e.g., [THUS]). One cannot 
rule out, however, that there exists a tricritical point at 
some finite length. 

In the original interpretation of [7| as a model with dis¬ 
tinct connectivity and dependency links, one can make 
either of them long ranged, or both. In |3S1|37||42] it is as¬ 
sumed that nodes are located on 2-d lattices with nearest 
neighbor connectivity links, while dependency links be¬ 
come long ranged. This is presumably not very realistic 
in most applications: Since dependency links are much 
more crucial than connectivity links, one would assume 
that they are the shorter ones in any realistic natural net¬ 
work. In the following we shall just make some comments 
on the model of [siiEaiig, and leave the more natural 
model of multiplex networks with long range connectivity 
links to a future publication. 
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2. Random Dependency Links 


As in [35], we first discuss the case where the depen¬ 
dency links are completely random (while, as we said, 
the connectivity links are the nearest neighbor bonds). 
The transition is obtained by varying the site occupancy 
q of one of the two dependency partners (i.e., any pair of 
dependency partners is present with probability q), and 
keeping all dependency links (i.e., p = 1). This is the 
r —>■ oo limit of the model discussed in the next sub¬ 
subsection. Following cascades, it was shown in |35| that 
the transition is first order in this limit, and closed for¬ 
mulas were given for the threshold qc and for the density 
of the viable cluster at threshold. 

As in the case of random locally tree-like graphs 
[M mi, these formulas are indeed obtained more eas¬ 
ily from consistency considerations without reference to 
any cascade dynamics. Let us denote by Sop{p) the or¬ 
der parameter (i.e. the density of the infinite cluster) in 
ordinary site percolation with site occupancy p (the sub¬ 
script “OP” stands for ‘ordinary percolation’). Let us 
now consider the percolation of one type of nodes, say of 
nodes of type A. Since dependencies are completely ran¬ 
dom, they are not correlated with the (non-trivial) order 
parameter fluctuations, and they just mean that a finite 
density of A nodes will be killed by not having a viable 
dependency partner. This fraction is exactly the same as 
the fraction of B nodes that are killed because they have 
no viable A partner. As a result, starting with occupancy 
p in ordinary site percolation leads to the same order pa¬ 
rameter as starting with occupancy q in the dependency 
model, provided that |3S] 

q=pxp/Sop{p) , (12) 


and Sop (p) is also the density of the viable cluster when 
starting with q. Plotting the r.h.s. of Eq. (121 versus 
p gives a convex function (see Fig. 25), the minimum of 
which is the critical value of q, 
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FIG. 25. (color online) Plot of p^/5'op(j3 ) versus p. The points 
are most precise near the minimum of the curve, where each 
point is obtained from « 10^ lattices of size 32768 x 32768. 
The statistical errors of these points are smaller than their 
sizes. 



gc = min pV'5'op(p) = 0.682892(5) (13) 

p 

with 

p* = argminp^/S'op(p) = 0.6418(2), (14) 

p 

and Sc = Sop{p*) = 0.6031(2) is the density of the giant 
viable cluster at threshold. The numerical values are in 
good agreement with the less precise estimates of [35] . 

The existence of a first order transition with these pa¬ 
rameter values was also confirmed by direct simulations. 
In Fig. 26 (main plot) we show the order parameter (the 
density of the largest cluster, averaged over all runs) plot¬ 
ted against q. Densities obtained by averaging only over 
those runs which did lead to a giant viable cluster are 
shown in the inset. As in the ER case, it was trivial 
to distinguish between runs with and without giant clus¬ 
ters, because all non-giant clusters had sizes 1 or 2 (ex¬ 
cept on very small lattices, where also clusters of size 3 


FIG. 26. (color online) Average density S of the largest viable 
cluster on 2-d lattices with random dependency links and de¬ 
pendency links between nearest neighbors, plotted against q. 
In the main figure, averages are taken over all runs (even they 
did not lead to a giant viable cluster, while the inset shows 
averages restricted to runs with giant viable clusters. 


were seen). Figure 26 confirms not only the values of qc 
and Sc but it also shows that the transition is hybrid, 
because the slope at threshold is infinite. The latter is 
a direct consequence of the fact that the minimum in 
Fig. 25 is quadratic. 

The dependence of the average SOS surface heights 
on q and on N = is also very similar to that for 
ER networks. As functions of q they show peaks that 
become narrower and higher with increasing L, and the 
peak values (more precise, the values at g = qc) show a 
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FIG. 27. (color online) Log-log plot of average SOS surface 
heights at q = qc for interdependent networks on square lat¬ 
tices with random dependency links. The straight line indi¬ 
cates a power law with exponent 0.556. 

power law 

{h{q,)) ^ (15) 

(see Fig. 27). This is in qualitative agreement with find¬ 
ings by D. Zhou (quoted in [12]), who obtained however 
an exponent 0.44, in clear disagreement with our value 
0.556. On the other hand, the scaling {h{qc)) 
is the same as that for ER networks (see Eq. (6)). 

3. Intermediate-Range Dependency Links 

Irrespective of all details, the results of the last sub¬ 
subsection show unambiguously that the transition is 
first order in the limit of infinitely long dependency links 
and nearest neighbor connectivity links, and the transi¬ 
tion occurs roughly at q < 0.7. This should be remem¬ 
bered if we now consider, following [311137], dependency 
links with intermediate range. 

In [m [37] the dependency links were random vectors 
X = {x,y) with ||x|| < r in the maximum norm, it was 
found that the model stays second order, as long as r < 8. 
At r = r* « 8 the transition becomes first order, and for 
r > 8 there are claimed to be two first order transitions 
with a metastable phase in between. The percolation is 
again implemented as site percolation, i.e. in one of the 
two lattices only a fraction q of nodes is open, while all 
connectivity bonds and all sites in the other lattice are 
open. 

We verified that there is an abrupt change of the tran¬ 
sition type at r = r*, although we found r* somewhat 
closer to 7 than to 8, see Fig. 28. The most conspicuous 
feature in the range 1 < r < r* is a dramatic increase 
of the surface heights. For all r the critical points are at 
those values of q where {h) changes fastest, which is to 



q 


FIG. 28. (color online) Average SOS surface heights plot¬ 
ted against q for interdependent networks on square lattices 
with dependency links of lengths < r, with r — 1 ,2, ... 8 . 
Each group of curves corresponds to the same value of r (in¬ 
creasing from left to right), while each curve within a group 
corresponds to a different L (with L = 512,1024,... 8192). 
For each r yf 7 the critical point is at the position where the 
curves for large L are steepest. Our estimate for r* is slightly 
larger than r* = 7. 


the left of the peaks for r < 7 and to its right for r = 8. 
For all r < 7 we found that the transition is not in the OP 
universality class, but in the universality class discussed 
in subsection III.B. This is most clearly seen for r = 1, 
where simulations are fastest (due to the smallness of {h)) 
and where also corrections to scaling are smallest. One 
might have anticipated that there are large scaling cor¬ 
rections for r = 1 because of the nearness to r = 0 where 
the model would be just ordinary site percolation, but 
this is not true. Indeed, we found critical exponents in 
perfect agreement with the estimates of subsection III.B, 
and with even somewhat smaller statistical errors. We 
do not show any plots as they are so similar to those in 
subsection III.B. 

For r > r* the transitions look very much like 
first order transitions, in agreement with the claims in 
[33I371I31]- If so, then the transition should be tricritical 
exactly at r *, but we were not able to measure tricritical 
exponents. As discussed in [SS] (ST] |42], what happens 
at r = r* is that fronts become unstable. Consider a 
scenario where open sites exist only in half space a: > 0, 
while all sites with x < 0 are closed. In ordinary perco¬ 
lation this would correspond just to a boundary where 
all clusters are cut off at x < 0. This case of percolation 
in the presence of boundaries has been studied in detail 
[saET]. In particular, the density of the giant cluster 
near the boundary is reduced and is governed by a new 
critical index, but otherwise the boundary has no effect 
on the bulk behavior. 

This is still true for r < r*, but not for r > r*. Al- 
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though we disagree with [33 |33 US] on details and on 
the theoretical treatment, we agree with them that nodes 
near the boundary are removed from the viable cluster by 
dependencies, which leads then to more removals, etc., so 
that such boundaries create moving fronts which finally 
destroy the entire viable cluster - unless q is large and r 
is sufficiently close to r*. If the front established in this 
latter case has a sharp profile and the density behind it 
is finite (which is clearly suggested by the simulations), 
one has a first order transition. 



q 


FIG. 29. (color online) Densities of the largest viable cluster 
for the model of [35] with r = 12, plotted against q for various 
system sizes. Here, averages are taken over all runs, including 
those with no giant viable cluster. 

This gives thus indeed rise to the “upper” of the two 
claimed transitions for r > r* (where qdr) > qcir*) ~ 
0.74) [3311371133] - with an important caveat discussed be¬ 
low. But we found no indication for the lower transition, 
where qc is supposed to be lower than the critical value 
at r*. In Fig. 29 we show the order parameter for r = 12 
(i.e., well above r*) for system sizes between L = 128 
and 8192. According to [33 [42], there should be a first 
order phase transition a.t q = 0.74. We see that there is 
a rather rapid cross-over at q = 0.74, if L « 2000 (which 
is presumably the size on which the claim of [33132] was 
based), but we can definitely rule out a real phase tran¬ 
sition within the range 0.68 < q < 0765. Indeed, extrap¬ 
olating to larger L we can be rather sure that no tran¬ 
sition occurs for any q < 0.8. Notice that qc = 0.753(1) 
for r = 8 (data not shown), so the transition point is 
definitely increasing with increasing r. We cannot say 
numerically whether it agrees with the transition defined 
via the propagation of fronts (is at g « 0.83 for r = 12), 
but this seems to be the only plausible option. 

Analogous results were obtained for r = 20 and r = 30 
(data not shown). Again we find rapid but nevertheless 
smooth crossovers at values of q which increase rapidly 
with L. In both cases the cross-over happens at the 
“transition points” seen in [33 132] when L « 10^, but we 


can rule out real phase transitions anywhere near these 
points. 



FIG. 30. (color online) Densities of the largest viable cluster 
for the model of |35| with r = 12 (same data as in Fig. 28), but 
conditioned on clusters which do have a giant viable cluster. 



q 


FIG. 31. (color online) Densities of giant viable clusters for 
r — 12, 20, 30, and oo. 

Further insight is obtained by plotting not the order 
parameter averaged over all runs, but by conditioning 
onto runs with a giant cluster. Again it is very easy to 
distinguish between runs with and without giant clusters 
(although the distinction is not as sharp as for r = oo). 
The results for r = 12 are plotted in Fig. 30. The per¬ 
fect data collapse shows that the difference between the 
curves for different L seen in Fig. 29 is entirely due to 
differences in the probabilities Pl with which a giant 
cluster is reached. The structures of the giant clusters 
themselves are completely independent of L. In Fig. 31, 
these data are plotted together with analogous results 
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for r = 20 and r = 30, and with the results for r = oo 
(obtained in the previous sub-subsection; see the inset in 
Fig. 26) extrapolated to L —oo. All four curves seem to 
become identical for large q, while they fan out system¬ 
atically for small q. This figure suggests that there would 
be very little dependence on r (except for very small g), 
if there were not a mechanism which would dramatically 
affect the probability for a giant viable cluster to occur. 



FIG. 32. (color online) Probabilities 1 — Poo that no giant 
viable cluster is formed, for r = 12. 

This mechanism is, as suggested in EH 02], the 
existence of voids in the initial configuration and their 
subsequent growth. This is also confirmed by plotting the 
probabilities that a giant viable cluster is not formed, see 
Fig. 32 for r = 12. For large q (i.e., if these probabilities 
are small), they increase with L as 

l-Poo^ c{q)L^, (16) 

suggesting that giant cluster formation is prevented by 
a rare extensive mechanism, i.e. by a mechanism whose 
probability to arise is oc This is of course true for the 
formation of voids in the initial configuration. We also 
see that these probabilities decrease, for fixed L, roughly 
exponentially with q. This is again what we would expect 
for the random formation of voids, with 

Pim)^—{l-qr (17) 

m 

being roughly the probability to form a void of m sites. 
Numerical estimates of void sizes based on this estimate 
and on Fig. 32 give void sizes which are typically half 
as large as those quotes in [35]. This is not unreason¬ 
able, given the fact that voids were assumed in |3S] to 
be circular disks, while they could indeed have different 
shapes. 

Thus we do agree with basic features proposed in 
[5511571 02] , but we do not agree with the existence of 
a second phase transition curve. We also would not call 


“metastable” the states for q slightly above this supposed 
transition curve, since metastability applies to systems 
subject to stochastic (not frozen) noise, and refers to 
states which are actually unstable on long time scales. In 
the present case we would rather speak of “conditional 
stability”, since the viable clusters are absolutely stable, 
but they arise only conditioned on particular initial con¬ 
figurations. 


4. A Final Remark on First Order Transitions and an 
Alternative Scenario 

As we said, the data suggest for r > r* a first order 
transition, and the location of this transition seems to 
be on the curve given in [55] [57102] . Although we have 
no direct numerical evidence to doubt this, there is in¬ 
direct evidence and theoretical arguments. As suggested 
in [55107102], the transition for r > r* is essentially a 
percolation transition, where the spreading of voids in¬ 
duced by dependencies becomes critical. In d > 3 there 
exist indeed tricritical points in generalized percolation 
models where the phase transition changes from second 
to hrst order Ellin]. There is however strong evidence 
[ElEE][59] that such tricritical points are absent in d = 2. 
This seems related to the Aizenman-Wehr theorem [50] 
on absence of first order phase transitions in random 2-d 
systems, although details are far from clear. Physically, 
it corresponds to the fact that phase boundaries cannot 
stay flat in 1-1-1 - dimensional isotropic random systems, 
but become crumpled on large scales. This is e.g. what is 
found in the T — Q random field Ising model [55] IM] ) and 
such interfaces seem always to be in the OP universality 
class, even if they look very smooth on small scales [5Tj . 

As suggested by Fig. 33, this is not true for simula¬ 
tions on lattices of size L_l x Ly with Ly ^ Lj_, and 
with zero initial densities in the regions x < X- = 20 
and a; > a;+ = Ly — 20. The latter enforce interfaces 
which would stay near x± ii q > qc, but move inwards 
for q < qc- We find, even for r very close to r*, that inter¬ 
faces dX q~ qc are very straight, with fluctuations much 
smaller than Lj_. For larger values of r (data not shown) 
the effect was much enhanced, and the interfaces were es¬ 
sentially straight lines. This is clearly a finite size effect, 
which is particularly enhanced by the specific choice of 
dependency links made in [55] 07] 02] : By choosing such 
links to be in a square region with sharp boundaries, fluc¬ 
tuations of the interface are strongly suppressed. 

We thus conclude that all simulations of 05107102] are 
far from the true asymptotic region. It might be that the 
positions of the critical points are nevertheless correct, 
but they may also be far off. In any case, the fact that the 
transitions for r > r* look like first order cannot be taken 
as evidence that they really are so. Explosive percolation 
[6] should be a warning that not all transitions that look 
first order also are so. In the present case, the length 
scale introduced by r should be enough to confuse the 
picture. 
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FIG. 33. (color online) Densities of giant viable clnsters on 
2-d lattices of size Lj_ x Ly = 2048 x 16384, with r = 9 and 
q = 0.772 (which is our best estimate for gc at r = 9). Shown 
are the 1-d projections of these densities for ten typical real¬ 
izations, multiplied with arbitrary factors to avoid overlaps. 
Due to local isotropy, we would expect that interface widths 
should become ~ L±_ in the large system limit. Actual inter¬ 
faces are much sharper. 


Indeed, we conjecture that the transition is in the uni¬ 
versality class of OP not for r < r* (as claimed in [35]), 
but for r > r* - where it just corresponds to ordinary 
percolation of voids. We have to admit that we do not 
have direct numerical evidence to support this claim, nor 
do we see any possibility to obtain evidence in favor or 
against it in the near future. 


IV. DISCUSSION 

The purpose of this paper was twofold: On the one 
hand, we proposed an efficient algorithm to simulate vi¬ 
able clusters in multiplex networks, and the percolation 
transitions related to them. This was motivated by the 
fact that exact results for this problem are obtainable 
only in mean field theory (i.e. on random locally loop¬ 
less graphs), and there exist ample indications that such 
predictions can be very misleading in real applications. 
Indeed, the algorithm presented in this paper is not only 
fast, but by mapping the problem onto a solid-on-solid 
model we were able also to discuss some non-trivial struc¬ 
tures which might become important for themselves. 

The other aim of the paper was to apply this algorithm 
in several simple networks and lattices, in order to test 
previous claims and to understand better the nature of 
the percolation transition(s) in this model. We found 
several surprises: 

• While the geometric structure of viable clusters on 
ER networks is indeed as predicted by mean field 


theory, it seems that the (pseudo-)dynamics of cas¬ 
cades is not. We speculated that this might be be¬ 
cause the clusters are supercritical during the cas¬ 
cades, but more work is needed to understand this. 

• On finite-dimensional lattices the model is not in 
the universality class of ordinary percolation, but 
represents its own new universality class. This is 
seen most dramatically in four and five dimensions, 
where OP would show continuous transitions with 
non-trivial anomalous exponents. Rather, we found 
a first order transition in d = 5, and a continuous 
transition with /3 = I in d = 4. Thus it seems as 
if the model has upper critical dimension du = 4, 
although the true mean field solution has a discon¬ 
tinuous transition. In d = 2 the results are less 
dramatic, but it seems clear that it is not in the 
OP universality class. In d = 3, finally, we can¬ 
not make a clear statement because of very strong 
corrections to scaling. 

• While the order parameter exponent is /3 = I for 
d = 4 as in other models at the upper critical di¬ 
mension, other exponents like the fractal cluster 
dimension and the correlation length are different. 
They seem to be simple rational numbers, but pre¬ 
cise estimates are difficult due to strong (logarith¬ 
mic?) corrections. 

• In general, all first order transitions are hybrid. 
Thus, the order parameter makes on ER networks 
and on 5-dimensional lattices not only a jump, but 
it also shows a power law. More surprisingly, in 
d = 5 there seems also to be a non-trivial scaling 
mass distribution of finite clusters. 

As regards geometric networks embedded in low di¬ 
mensions with semi-local links, we partly confirmed the 
scenario found in [35]. We do not find the double phase 
transitions claimed there when the lengths of these links 
are are above a (tri-)critical value. We do find that there 
exists such a (tri-)critical value above which the transi¬ 
tion seems to be first order, but we give strong evidence 
that this is related to very large finite size corrections. 
More importantly, we claim that the model studied in 
[351I3Z1 I42j is very unrealistic in assuming dependency 
links to be much longer than connectivity links. Whether 
a model where this relation is reversed shows similar be¬ 
havior is an open question. 

Other open questions concern the behavior of multi¬ 
plex networks with > 2 types of links. We have not yet 
tried to extend our algorithm to this case. Even less 
clear is the behavior in case of non-mutual (asymmetric) 
dependencies, which cannot be mapped onto multiplex 
networks at all. The same is true for directed multiplex 
networks. Finally, the true asymptotic behavior in the 
2 -d case with long range dependency links could, eventu¬ 
ally, only be solved analytically or by means of a different 
model realization which avoids the introduction of a new 
length scale. 
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